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PREFACE 

This  report  summarizes  Pacific-Sierra  Research  Corporation's 
(PSR)  work  in  developing  a  practical  numerical  method  for  computing 
radio  field  strengths  in  regions  of  strong  focusing,  using  input 
from  a  standard  ray  tracing  program.  The  computer  program  was 
designed  for  use  with  the  Jones-Stephenson  high-frequency  ray 
trace  code. 


FKEGEOUO  tiOK  MJJK-MOT  HliCD 


I .  INTRODUCTION 


This  report  presents  a  method  for  using  ray  tracing  to  calculate 
radio  field  strengths  in  strong-focusing  regions  of  the  ionosphere 
or  troposphere.  Ray  tracing  is  an  established  technique  for  predict¬ 
ing  radio  field  strengths  in  regions  without  strong  focusing.  Because 
of  the  versatility  and  numerical  efficiency  of  ray  tracing,  it  is  a 
desirable  technique  for  computing  radio  fields  in  strong-focusing 
regions,  using  a  state-of-the-art  ray  tracing  program  such  as  that 
of  Jones  and  Stephenson  [1975].  This  report  presents  a  practical 
numerical  method  for  doing  so. 

The  procedure  for  computing  radio  fields  in  stratified  media  from 
ray  tracing  outside  caustic  regions  uses  the  well-known  correspondence — 
developed  by  Booker  [1939],  Budden  [1961a],  and  others — between  ray 
trajectories  and  the  first-order  asymptotic  approximation  to  the 
angular  spectral  representation  of  the  field  components.  Budden  [1976] 
and  others  show  how  the  phase  and  amplitude  of  fields  are  directly 
related  to  the  path-integrated  refractive-index  variation  and  the  ray 
curvature,  both  easily  computed  with  a  ray  tracing  program  in  regions 
where  neighboring  rays  do  not  cross.  That  correspondence  solves  the 
radio  field  problem  for  stratified  media  outside  caustic  areas. 

Breakdown  of  the  first-order  stationary-phase  result  near  caustics 
has  led  to  the  development  of  higher  order  asymptotic  formulas  that 
uniformly  interpolate  the  field  through  caustic  regions  and  can  extend 
field  strength  calculations  based  on  ray  density  into  strong  focusing 
regions  [Maslin,  1976i;  Budden,  1976].  Unfortunately,  asymptotic 


methods  have  two  drawbacks  that  make  them  inconvenient  for  use  with  a 


ray  tracing  program; 

1.  Different  closed-form  expressions  involving  special  functions 
are  needed  for  different  degrees  of  focusing — e.g.,  Airy  func¬ 
tions  for  caustics  and  Pcarepy'c  [1946]  function  for  cusps. 
Further,  only  the  Airy  functions  are  both  well  known  and 
convenient  for  numerical  calculation. 

2.  The  contributions  from  all  rays  intercepting  each  field  point 
must  be  calculated — a  numerically  cumbersome  procedure  re¬ 
quiring  a  search  to  identify  the  specific  plane  wave  compo¬ 
nents  intercepting  a  given  field  point. 

Because  of  these  limitations,  only  simple  refractive-index  profiles 
can  be  studied,  wasting  much  of  the  potential  power  of  the  ray  trace. 

An  alternative  to  asymptotic  methods  is  to  evaluate  the  integral 
expressions  for  the  field  numerically  using  phase  information  derived 
from  the  ray  trace.  Thus,  a  broader  range  of  refractive-index  pro¬ 
files  could  be  analyzed  by  a  single  method;  cusps  and  even  higher 
order  focusing  effects  would  require  no  special  analysis.  However, 
this  method  has  two  potential  difficulties: 

1.  The  Integrands  are  highly  oscillatory  and  standard  numerical 
integration  techniques  fail. 

2.  For  points  in  the  vicinity  of  caustics,  more  accurate  phase 
differences  between  neighboring  plane-wave  spectral  compo¬ 
nents  are  needed  than  the  ray  trace  can  supply  directly 
through  its  path  integration  of  the  refractive-index  varia¬ 
tion.  That  limitation  is  simply  due  to  a  loss  of  numerical 
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precision  and  disappears  if  the  ray  directions  through  a 
given  field  point  are  widely  separated — as  indeed  they  are 
when  the  first-order  calculation  is  valid. 

The  numerical  method  presented  here  overcomes  tliese  difficulties. 

Section  II  briefly  sunmarizes  the  plane-wave  angular  spectral 
representation  of  radio  fields,  including  Kuddrn':^  |197f>)  generaliza¬ 
tion  to  allow  the  phase  integral  method  to  ne  applied  uniformly  through¬ 
out  reflection  ragions.  Section  III  discusses  a  method  for  using  ray 
trace  data  to  accurately  estimate  the  phase  differences  between  neigh¬ 
boring  plane-wave  components.  Section  IV  shows  that  the  highly  oscil¬ 
latory  spectral  integrations  can  be  performed  by  a  numerical  quadrature 
technique  developed  independently  by  Woodie  [1976]  and  Ik2i‘akat  [1976]. 
Field  calculations  using  our  method  are  compared  with  asymptotic  re¬ 
sults  of  Maslin  [1976b]. 

The  only  potentially  important  feature  lost  with  our  approacii  is 
the  ability  to  include  eva'nescent  wave  components  in  the  spectral 
integral;  they  must  be  excluded  if  the  ray  trace  can  propagate  only 
real  rays.  That  limitation  precludes  treating  problems  such  as  leak¬ 
age  through  layers  at  frequencies  slightly  below  the  penetration  fre¬ 
quency.  Section  V  gives  a  partial  solution  to  the  leakage  problem, 
generalizing  Buddan’i:  single-reflection  height  interpolation  to  the 
case  of  two  arbitrarily  close  roots  of  the  /^ookwy  q  function. 

For  convenience,  mathematical  details  of  the  phase  determination 
and  numerical  integration  procedures  are  given  in  two  appendixes. 
Throughout  our  analysis,  we  assume  an  isotropic  propagation  medium 
varying  in  only  one  dimension — wlilch,  however,  can  he  either  height 
(flat  earth)  or  radius  (spherical  earth). 


II.  UNIFORM  AN('.lll.AR  SPECTRAL  REPRE.SF.NTATION 


Tlu'  basis  of  our  analysis  is  Hudeicn'r  [197f)l  integral  expression 
lor  a  Field  eomponent  F  that  is  uniformly  valid  tlirougbout  caustic  and 
reFleclion  regions: 


written  a.s  a  superposition,  weighted  by  the  function  G,  of  plane 
waves  with  complex  direction  cosines  ,  S^,  C  in  free  space.  We 
use  Cartesian  coordinates  x,  y,  z,  with  z  denoting  altitude;  k  is 
the  free-space  wave  number.  The  function  q  is  defined  as 

q^  =  n^(z)  -  sj  -  .S^  ,  (2) 

where  the  complex  refractive  index  n  Is  t.aken  to  lie  a  function  of 
lieight  z  only.  Following  f,  is  defined  as 


Ai  is  Airy's  integr.il,  .ind  x^(.S|,  .S^)  is  the  reflection  height 

for  the  pl;ine  wave  component  labeled  S^,  f'uddrn 


(19761  provides  the  proper  Integration  contour  and  phase  choices 


that  make  (1)  unambiguous.  When  »  1,  tlie  factor 


Aifrj  exp  ( -ik  I  clzq 


/ 


(4) 


assumes  the  asymptotic  form 


7.  \  /  y. 

exp  I  -ik  I  dzq  )  +  i  ‘•'xp  (  - 

•'0  /  \  "O 


y-  \ 

ik  <J/.c| 

j 


(■i) 


representing  a  sum  of  upgoing  and  downgoing  waves.  Tlu*  |)aLli  integral 
over  q  in  the  second  term  is  understood  to  eircU*  tlie  rc'f  I  I'ct  ion 
height  2q  clockwise  in  the  analytical  continuation  of  heiglit;  on  tlie 
downgoing  path,  q  is  chosen  negative.  Although  it  is  valid  to  re¬ 
place  the  Airy  function  factor  (4)  with  tlie  asympt  ol  i  i-  limit  (5)  even 
at  caustics  when  |5|  1  1976a],  the  Airy  form  is  capable 

of  including  plane  wave  components  near  a  reflect itin  level,  where 
q(Sj^,  S2)  vanishes.  for  completeness  we  prefi>r  to  retain  the  general 
form  (1)  although,  for  most  spectral  components  appisiring  in  the  in¬ 
tegral,  the  asymptotic  approximation  using  Ci)  is  both  adequate  and 
faster  to  compute.  The  asymptotic  form  is  in  fact  the  basis  for 
the  higher  order  as}anptotic  approximat Ions  of  fields  near  caustics 
and  cusps  constructed  by  MaP.Vi.n  l\97(>b].  His  results  contain  Airy 
functions  only  because  ot  his  use  of  the  uniform  asymptotic  method 
of  Chester  eb  al.  [1957]. 

The  spectral  weighting  function  C  is  determined  by  both  the 
field  component  represented  by  F  and  the  transmitter.  For  most  appli¬ 
cations,  the  angular  spectral  dependence  of  C  will  be  fairly  simple. 
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1 1‘,  llic  spi'fliMl  weight  in>>  of  tlio  y-component  of  tho  olcc- 
Lric  field  in  tlic*  magnoL  ic  dipole  example  studied  by  Mac, Lin  11976/)] 
and  Hnddfn  [19761  is  represented  by 

C(Sj.  S^)  =  S^/C  .  (6) 

It  is  thus  eonvenii’iU  to  expand  G  in  Ijarmonies  of  its  azimuthal  dc- 
pendenne  and  treat  each  component  separately.  That  is,  we  expand  G 
as 


0f> 


(7) 


wliere 


2ir 


«n(^)  =  2ir 


d(})  G(S 


1 


S2) 


-in(]> 


(8) 


and  =  S  cos  ({),  =  S  sin  (f).  Substituting  (7)  into  (1)  and  inte¬ 

grating  over  ({)  produces 


l•(x,  y,  z)  =  2ir 


ini]) 


n=-o 


(I 


dS  S  g^fSXC/q)'"^^  A(S,  z)J^(kSp) 


(9) 


where  x  =  p  cos  >]<,  y  =  p  sin  \j),  and  A(S,  z)  denotes  either  the  Airy 
function  factor  (4)  or  its  asymptotic  approximation  (6).  For  points 
not  directly  over  tlie  transmitter,  the  Bessel  functions  may  be  approx¬ 


imated  as 
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“(jiss)  ”  I  - 1 )  •  <’“> 

Specializing  to  the  magnetic  dipole  represented  by  (6)  for  the  plane 
y  =  0,  using  (10)  and  the  asymptotic  form  of  A,  yieltls 


F(x.  0,  z) 


(2'rr/kx)^^^  y*  dS  (Cq)"’^"^ 


(H) 


for  the  upgolng  wave.  An  analogous  expression  results  for  tlie  re¬ 
flected  wave.  Equation  (11)  here  is  equivalent  to  Eq.  (7)  of  Mnclin 
[19762?]  and  Eq.  (23)  of  HiAdden  [1976). 
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COMl’DTA'nON  OK  IMIASK.  KKOM  KAY  'I'KACM'  DATA 


I'or  siif  f  icicnl  ly  .simple  re  f  rae  t  ivt— index  profiles,  the  phase 
integrals  such  ns 

V. 

dzq(Sj^,  S^,  z) 


needed  in  (9)  are  most  simply  computed  with  nnmc'rical  or  analytic 
methods.  For  complicated  media,  however,  .such  evaluation  becomes 
tedious — especially  when  ducting  with  multipU'  turning  point.s  is  in¬ 
volved.  On  the  other  hand,  ray  tracing  techniques  have  been  found 
useful  In  producing  field  estimates  through  ray  density  calculations 
for  quite  general  ionospheric  or  tropospher !<•  models  [Wnnij,  1958; 
l<ii(l(ii;n  tmd  Tcffii ,  197]].  The  principal  d  is.idvantage  of  ray  density 
calculations  is  their  breakdown  in  caustic  regions.  Nonetheless,  the 
identification  of  rays  with  spectral  w.ave  components  is  the  basis  for 
;i  simple  moans  ol'  using  the  ray  trace  to  estimate*  the  phase  inte¬ 
grals  required  in  the  HuddrnlMai'lin  theory  discussed  above. 

The  connection  between  ray  tracing  and  phase  integration  can  be 
established  through  the  Hamiltonian  formalism.  The  following  are  the 
canonical  equations  [llar,(:L(jrni\’  and  lUir.al  <jnov<' ,  19f>0)  for  the  ray 
trajectory  spi’cified  by  position  x,  y,  z  and  direction  cosines  S^, 

S^,  q,  where  .S^  +  +  q  =  n  (/.) : 


d_x  ^  llll  iX  =  M 

dT  il.Sj  ’  dT  ~  ’  ilT  °  dq 


...Ay 


(12) 
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and 

^1  _  _  ^  ^  M  iia.  =  i'* 

dT  “  3x  ’  dr  “  Dy  ’  dt  "  Dx  ' 


Here  T  parameterizes  ttie  ray  trajee.tory  and  ean  be  taken  to  be  the 
arc  length  along  the  ray.  The  Hamiltonian  H  can  be  specified  in  a 
number  of  ways;  a  convenient  choice  is 


H 


2  2 
+  S2  +  q 


0  . 


(14) 


The  ray  trajectory  results  from  substituting  (14)  into  (12)  and  (13), 
then  solving  them  either  numerically  or  analytically.  Sinc-e  .S^  and 
S2  are  conserved  along  the  ray  trajectory  (Snell's  law),  the  spectral 
component  labeled  S^,  S^  can  be  identified  with  a  unique  ray.  From 
elementary  calculus 


dx  ^  ^  ^  ^  ■  dy  ^2 

dz  “  dT  /  dT  "  as^  /  Dq  "  q  dz  "  q' 


and  integration  gives 


(15) 


(Ifvt) 


(16/0 


where  (14)  was  used  for  the  last  step.  Kquat  ions  (16)  are  identical 
to  Eq.  (5)  of  Mar.l.iri  (1976/;1  and  are  the  ray  trai'ing  equations  of 
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Haokcr  [1939].  Manlin  nnd  Itiuklrn  derive  them  by  applying  the  first- 
order  stationnry-phase  method  to  the  angular  spectral  representation 
(1) .  The  important  point  is  that  (16)  remain  valid  when  the  ray 
passes  through  a  caustic,  even  though  the  associated  first-order 
stationary-phase  approximation  to  the  field  breaks  down.  Aside  from 
a  constant,  the  phase  integral 

dzq(Sj,  S^,  z) 

is  recovered  by  integrating  (16)  with  respect  to  or  S^.  The  inte¬ 
grands  X  or  y  are  produced  using  ray  intercepts  on  the  plane  at  height 
z  as  a  function  of  Sj^,  S2.  Appendix  A  details  this  construction. 

Because  this  method  of  constructing  the  phase  from  the  ray  trace 
uses  the  ray  intercepts  and  not  the  path-integrated  refractive-index 
variation  directly,  it  provides  a  much  more  reliable  estimate  of  the 
small  phase  differences  between  neighboring  plane  wave  components. 

This  is  essential  to  the  success  of  the  approach. 

Although  (J)  is  in  general  complex,  only  its  real  component  can  be 
recovered  by  this  method,  which  employs  a  real-valued  ray  tracing 
program  sucli  as  that  of  Jcmi  n/.'Uci’hfru'cn  [1975].  Tims,  a  primary 
assumption  of  our  model  is  that  evanc.'scent  wavi's  rontrihute  negligibly 
to  the  field  represeiUetl  by  (9).  Nevertheless,  weak  collisions  may 
be  ineliided  In  this  approach  by  formally  tracing  the  rays  as  if  they 
were  real,  then  adding,  the  p.'ith-iiUegrated  attenuation  computed  by 
the  ray  trace  to  the  spectral  weighting  function  for  each  wave  com¬ 
ponent  [/iooknr,  1939). 


({)(S| ,  S^,  7.)  = 


I 
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The  neglect  of  evanescent  waves  appears  to  be  valid  as  long  as 
the  spectral  weighting  function  G  is  continuous.  For  G  continuous, 
the  absence  in  the  integral  of  spectral  components  with  a  slowly  vary¬ 
ing  phase  in  the  shadow  region  beyond  the  caustic  leaves  only  rapidly 
varying  components,  leading  to  destructive  interference  and  a  decaying 
field.  For  an  infinitely  distant  transmitter,  however,  there  is  only 
a  single  plane-wave  component  S^;  G  tlius  becomes  a  delta  function 
6(S  -  Sq)  and  our  model  would  not  provide  information  about  the  field 
above  the  reflection  height  of  the  plane  wave  component  S^. 

Because  q  is  double-valued,  (16)  can  he  taken  to  rc’presc*nt  Intel — 
cepts  of  either  the  upgoing  or  the  reflected  parts  of  the  ray  tra¬ 
jectory.  In  the  Latter  case,  the  integral  is  understood  to  run  up 
to  Zq,  then  to  return  to  z  with  q  changing  sign  as  it  passi-s  through 
zero.  Thus  we  d<>fine  the  ph.ise  of  the  up-  and  downgoing  rays  and 
(^d  as 


^2’ 


7,)  = 


Z 

f. 


dzq(S^,  S^,  z) 


(17) 


and 


♦d<"] 


'2’ 


dzq(Sj , 


7.) 


2’ 


.  • 


dzq  . 


(18) 
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Tho  ray  trnjortory  is  symmc'trical  about  its  n-f  J  oil  ion  point,  so  (18) 
can  bo  written 


2’ 


/.)  +  2 


/ 


dz|q| 


(19) 


Intogration  of  tlio  up-  and  downgoinp,  ray  intoroopts  produces  (J)^^  and 
4’j  on  tlio  plane  at  z  only  up  to  i-onstants;  tbe  proper  phase  differ- 
iMice  between  the  up-  and  downgoin};  waves  is  produced  by  requiring  the 
phases  to  agree  at  the  reflection  point. 

Suhst 1  tut Ing  those  results  into  (3),  we  obtain 


■  [i  ""♦u  -  ♦„>] 


2/3 


(20) 


and  (9)  becomes 


l'(x,  y,  z)  =  2it 


^  c'"*  f\ 

- -  n 


dS  S  g  (S)(C/q)'^^-  2f/^^  Ai(f,) 

n 


exp  l^-l  ^  (<^J  +  <}>^^)j  .l^(kSp)  . 


(21) 


Using  (10)  for  points  nut  directly  over  the  transmitter  yields 


l’(x,  y,  z)  =.e‘’'^‘  (2ll/kp)'^^ 


ns-on  •'-1 


.  .  iff 


1/2  ,1/A  (  ..  r.,  .  1 

r  A  I  /  r  %  |»vn  ' 


X  2fl‘'^  f/'"'  Ai(r.)  exp  I  2  ^“^d  ^  | 


^d  v] 


(22) 


In  equations  (21)  and  (22),  the  actual  liUe>;rat  ion  limits  are  deter¬ 


mined  by  tile  rays  tliat  inter<-epl  the  plane  al  ■/. .  l•■or  |  T.  |  ••  1,  the 

factor 

271^^^  exp  ((j)j  +  (23) 

becomes 


exp  (-ik(j)^^)  +  i  exp  (-ik(j)^)  .  (24) 

in  practice,  Huddcn  and  ftn'ty  [1971]  show  ]f^|  need  only  bc“  . 3  to 
make  (24)  a  good  estimate  of  (23).  Following  Hudden  [1961],  ^  is 
computed  as 

.  2/3 

^  (  k(4,d  -  4>„)  .  (25) 

and 


^1/4  irl'/^ 

=  !C!  f 


(26) 


These  are  the  proper  phase  choices  for  producing  the  asymptotic  form 
(24)  from  (23)  well  below  the  reflection  height. 
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IV.  NUMKKJCAI. 

Aithoiigli  our  mi'Lhod  i.s  suitable  for  .'i  wide  range  of  rcalistie 
profile.^,  we  iMiistrntc  it  by  appl  ic.'ition  to  two  relatively  simple 
ionospheric  profiles  studied  by  Mat'fin  |l976,;,/i|,  allowing  our 
numerical  results  to  he  compared  with  his  asymptotic  results. 

I.INKAR  rONOSPHKRIC  PROKILK 

Rirst  consider  the  linear  ref r-ictlve- index  height  profile  used 
by  Mdi'  /  in  [  1 976<f ,/)  1  : 

n  =  I  ,  z  <  h  , 

=  I  -  oi(z  -  h)  ,  z  >  h  ,  (27) 

where  a  and  b  are  I'onstant.s.  Hecause  of  the  simplicity  of  Die  model, 

an.'ilytic  ex))ress  iiins  for  tlu*  ray  intercepts  on  an  arbitrary  horizontal 
pl.'ine  at  height  z  ;ire  easily  found  from  (2)  and  (lb);  in  the  plane 

y  =  0, 

x(S)  =  zS/C  ,  z  <  h  Iipgoing  , 

=  (2li  +  4(;^/ot  -  z).S/C  ,  z  <  h  downgoing  , 

=  liS/C  +  -  ni?.  -  h)  J  2.S/a  ,  z  ^  li  upgoing  , 

=  hS/(!  +  ^(1  +  ~  -  li)  ^  2.S/(r  ,  z  ■  b  downgoing  . 

(28) 

l•’lgure  1  plots  x(S)  as  a  function  of  .S  for  various  heiglits  z,  using 
h  =  100  km,  a  =  0.002  km  *.  Tor  graphical  clarity,  the  intercepts 


Ray  intercept,  x  (km) 


Fig,  1--Ray  intercepts  on  horizontal  planes  at  various  heights  z 
the  linear  ionospheric  model  (27),  with  h  =  100  km, 
a  =  0.002  km-1.  Caustic  loci  are  shown  as 
dashed  lines  that  meet  at  two  cusps. 


have  been  left.  Incomplete  for  z  =  110  and  120  km.  Dashed  lines  plot  the 
loci  of  the  caustics.  The  loci  meet  at  two  cusps — one  below  the  ground 


at  X  ~  1360  km,  tlie  other  at  a  height  of  abt)iit  11 '3  km  at  x  ~  800  km. 

A  comparison  of  our  figure  with  Figure*  3  of  |1976al  helps 

to  I'larify  the  conf  igur.'ttion  of  the  ray  paths.  For  each  lu'ight 

7.  in  our  figure,  the  up-  and  downgoing  rays  form  the  lower  and  upper 

br.nnches  of  a  closed  loop,  respectively.  The  branches  connect  at  the 

reflection  vnliu*  S  (/.);  rays  with  S  >  S  (z)  do  not  reach  the  height  z. 
I'  c 

The  figure  also  shows  that  the  caustics  are  associated  only  with 
the  downgoing  branches. 

For  the  magnetic  dipole  repre.sented  I  v  (6)  (inly  the  n  =  ■•  J  terms 
contribute  to  (22).  In  the  plane  y  =  0, 


F(x,  0,  z) 


ill /A 

(.1 


1  /'> 

(2ti/kx) 


exp  I  - I k  Sx 


+  2 


♦  v] 


(29) 


with  r,  given  by  (2')).  Here  F  represents  the  y-component  of  the  elec¬ 
tric  field.  The  phase  integrals  (})  .  and  (b  can  be  found  either  an- 
alytically  (in  this  case)  or  through  the  ray  tracing  proceiluri*  dis¬ 
cussed  in  Sec.  III.  To  conifKirc  with  M.;rl  in'.-  |1976,'  )  a.symp  t  ot  ic  cal¬ 
culation,  we  set  z  =  0,  k  =  Ati  km  Owing  to  the  very  large  phas(' 
oscillations  in  (29),  its  numerical  evaluation  requires  special  cari’. 
Appendix  B  describes  a  numerical  quadrature  method  developed  inde¬ 
pendently  by  Woodic  [1976]  and  l^anakat  [1976],  the  piecewise  linear 
phase  (PLP)  algorithm.  Although  more  standard  methods  such  as 
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the  fast  Fourier  transform  (FFT)  could  be  as  efficient  as  tlie  PLP 
algorithm,  we  liavi'  ohtainejl  consistently  good  i  i-su  1 1  s  using  tlii’  l.ittiT 
with  highly  oscillatory  integrals.  Its  application  lu-re  with  100  suh- 
intcrvals  produced  the  result  slu>wn  in  Figuri-  2  tor  tlu-  modulus  of  the 
field.  Our  figure  agrees  clost'ly  with  Figure  1  of  hUiy.lin  [I97h/<1, 
even  in  the  vicinity  of  the  i-.austics  .-it  1189. ft  ami  122b.  7  km. 

HYPERBOLIC  SECANT  PROFILE  W^I TH  (p U. I S  10{« 

A  well-known  model  of  an  ionospheric  layer  is  the  hyperbolii- 
secant  profile 

2 

_  j  _  .sech^  la(z  -  1  ,  (30) 


where  a,  a,  z  are  constants.  For  Z  =  0,  the  ray  intercepts  in 
m 

the  plane  y  =  0  are  given  by  Eq.  (32)  of  hUu'.lin  [197b(t): 


1  2  2  2 ' 

„  )C  cosh  {-0.7.  )  -  a  I  -  C  sinh  (-txz  ) 

S  -  _  _ _  m _ _  _) _  _  _  _ m 

i  “2  "*'2"  . . 2!''^^  ' 

+  !c  co.sh  la(z  -/.))-  a  [  -  C  sinh  Ia(/. 

L  I  ) 


Figure  3  plots  x(S)  for  various  z,  using  a  =  0.9,  u  =  0.05  km 

7,  =  100  km,  and  ignoring  collisions  (Z  =  0).  For  this  cIkmci’  of  t  lu' 

m 

parameter  a,  r.iys  having  S  <  =  0.A3  penetrate  tlu-  layer  and  an- 

not  reflected  at  ;iny  lu-ight.  The  interpretation  of  Figure  3  is  similar 
to  that  of  Figure  1;  the  caustics  are  sitown  as  dashed  lines  meeting  ;it 
a  single  cusp  ;it  z  91  km.  Again,  for  graphical  clarity,  two  of  the 
iipgoing  branches  have  been  left  incomplete. 


--Range  dependence  of  the  electric  field  modulus  at  the  ground  for  the 
linear  ionospheric  model  (27),  with  h  =  100  km,  a  =  0.002  km"!. 


To  rest  Clio  phase- integral  cstlination  teclinique  described  in 
Sec.  Ill,  we  used  I  lie  Joniu’./HLephmaon  (19751  ray  tracing  program  to 
compute  ray  intercepts  x(S)  on  the  horizontal  plane  at  height  88  km 
for  both  up-  and  downgoing  rays  with  S  values  between  0.47  and  0.65. 
The  fits  were  adjusted  to  yield  <j)^  =  for  .S  =  0.65,  corresponding 
to  the  ray  grazing  the  plane  at  z  =  88  km. 

Collisions  were  modeled  using  the  procedure  suggested  by 
Hookaf  (1939]  of  propagating  the  rays  as  though  they  were  real,  but 
modifying  the  spectral  weighting  function  G  by  the  path-integrated 
aCConuation  for  each  ray  traced.  This  was  easily  done,  since  the 
Jonci'.  It'll  <  phcni'.nn  program  can  compute  the  attenuation.  For  comparison 
with  Mdr.lin  (1976/;],  constant  colli.sions  given  by  the  parameter  Z  of 
0,  0.0005,  0.0018,  and  0.0045  were  used.  Application  of  the  PLP 
algorithm  to  the  magnetic  dipole  in  the  plane  y  =  0  produced  the 
curves  shown  in  Figure  4  for  the  modulus  of  tlu’  y-component  of  the 
field  as  a  function  of  range.  Oiir  figure  agrees  well  with  Figure  6 
of  Mar, Lin  [1976/»1.  The  two  figures  show  the  same  reduction  in  oscil¬ 
lations  for  larger  collision  frequencies. 


-Range  dependence  of  the  electric  field  modulus  at  a  height  z  of  88  km  for  the 
linear  ionospheric  model  n2  =  1  -  a2  sech2  [a(z  -  z^,)].  with  a  =  0.9, 
a  =  0.05  km~l ,  Zm  =  100  kra.  The  different  curves  represent 
constant  collision  frequencies  specified  by  the  parameter  z- 
A,  0;  8,  0.0005;  C,  0.0018;  D,  0.0045. 
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V.  UNIFORM  ANCULAR  SPECTRAL  REPRESENTATION  FOR  AN  EVANESCENT  LAYER 


The  foregoing  treatment  of  fields  in  regions  of  strong  focusing 
using  ray  tracing  techniques  was  based  on  Ihuiden's  [1976]  approximate 
solution  (1)  to  the  wave  equation,  uniformly  valid  in  the  neighbor¬ 
hood  of  a  single  reflection  height.  His  solution  is  a  good  approxi¬ 
mation  to  the  exact  wave-equation  solution  throughout  the  reflection 
region  if,  first,  higher  order  derivatives  of  the  refractive  index 
with  respect  to  altitude  are  well  behaved;  and,  second,  the  reflection 
height  is  sufficiently  remote  from  other  roots  of  Hooker' 's  [1939] 
quartic  function  q.  The  second  condition,  however,  excludes  problems 
such  as  the  leakage  of  tropospheric  ducts. 

Considerable  work  has  been  done  in  modeling  duct  propagation 
using  the  waveguide  mode  approach  [Hook,  1950;  Fook,  Weinstein,  and 
Belkina,  1956,  1958;  Budden,  1961/>;  Wait,  1962;  Chang,  1971;  Papper't 
and  Goodhavt,  1977,  1979].  However,  no  one  has  generalized  the 
angular  spectral  representation  to  include  the  case  of  partial  trans¬ 
mission  through  a  thin  evanescent  layer  in  which  there  are  two  neigh¬ 
boring  roots  of  q.  That  extension  is  needed  to  allow  our  analysis 
to  be  applied  to  superrefraction  in  leaking  tropospheric  ducts.  We 
outline  the  construction  of  the  generalized  spectral  integral  holding 
uniformly  through  evanescent  regions  and  show  that  asymptotic  limits 
of  the  expression  produce  the  correct  results  for  thick  layers. 

It  is  well  known  that  any  field  component  F  in  free  space  can 
be  written  as  a  superposition  of  plane  waves  as 


T 
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F(x,  y,  z)  = 


If 


dS,  dS„  G(S 


1’ 


S2) 


exp 


[- 


lk(S^x 


^2>' 


+  Cz) 


(32) 


In  a  medium  varying  in  z  only,  the  x-  and  y-dependence  of  a 
given  spectral  component  of  F  in  (32)  is  unchanged  but  the  z-dependence 
is  determined  approximately  by  the  equation 

3^F  ^  ,  2  2„  _ 

— X  +  k  q  F  -  0  .  (33) 

3z 

Additional  terms  involving  derivatives  of  n  may  be  present  in  (33), 
depending  on  the  choice  of  transmitter  and  field  component;  they  are 
neglected  here.  We  wish  to  solve  (33)  approximately  for  the  case  in 
which  there  are  two  arbitrarily  close  roots  ^  of  q.  For  con¬ 
creteness,  we  assume 


2 

q  >  0  ,  z  >  z^,  z  <  Zj^  , 
2 

q  <  0  ,  z^  <  z  <  Z2  • 


(34) 


This  choice  corresponds  to  the  case  of  an  evanescent  region  between 
Zj^  <  z  <  Z2  with  real,  propagating  rays  outside.  The  complementary 
problem  of  taking 


2 

q  >  0  ,  <  z  <  Z2  , 

2 

q  <  0  ,  z  >  Z2,  z  >  , 


(34a) 


represents  an  elevated  duct.  Because  the  latter  case  has  been  well 
researched  [Fock,  Weinstein,  and  Belkina,  1956,  1958;  Wail,  1962], 
we  concentrate  on  the  evanescent-layer  problem. 


For  large  k,  (33)  is  a  special  case  of  a  more  general  equation 
studied  by  Lan/jer  [1959J.  The  essence  of  his  method  is  to  change 
variables  to  produce  an  equation  resembling  Weber's  [1869]  equation 
asymptotically  as  k  The  change  of  variables  is  defined,  first, 

by  the  equation 


with 


Q  2 


and,  second,  by 


(35) 


(36) 


(^2  _ 

F(z)  =  .  (37) 


in  terms  of  which  (33)  becomes 


(38) 


Here 


7. 


I 


^  _  (d  -  1)^^^  2Q 
dC  q  H  ’ 


(39) 


with  analogous  expressions  for  z"  and  z"’ .  Following  hanger,  the 

2  1/2 

arguments  of  q  and  (^  -  1)  in  (35)  are  defined  as 
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0,  z  <  and  5  ~1  . 

7r/2,  Zj  <  z  <  Z2  and  -1  <  ^  <  1  , 

0,  z  >  Z2  and  5^1* 

(AO) 

Equation  (35)  thus  can  be  used  to  define  a  mapping  between  heights 
throughout  the  region  of  the  layer  and  the  real  line.  For  sufficiently 
high  frequencies,  the  last  two  terms  in  (38)  are  small  relative  to  the 
first  two;  dropping  them  and  using  the  substitution 


^  ^  /■t-2  ,.1/2 

arg  q(z)  =  arg  (^  -  1) 


t 


W(t)  =  uiO  . 


produces  the  following  form  of  Weber's  equation: 


(Al) 


Besides  Weber  [1869],  many  others  have  studied  solutions  to  this  equa¬ 
tion,  including  Damoin  [19A9],  Miller  [1952],  Olvcr  [1959],  and  Erdehji, 
Kennedy,  and  MaGreyov  [195A].  Abrarnowit::  and  Slegun  [1965]  summarize 
those  previous  analyses.  Using  the  normalization  and  notation  in 
Abramowitz  and  Stegun,  two  independent  solutions  of  (A2)  are  given 
by 

Physical  requirements  dictate  the  proper  combination  of  these 
solutions.  For  kQ/ir  >>  1,  we  require  the  result  to  agree  with  Huddcn’s 
[1976]  uniform  asymptotic  result  for  z  near  z^.  For  z  >>  7.2,  the 
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solution  should  represent  an  outgoing  plane  wave.  As  shown  below, 
those  limits  can  be  achieved  by  choosing  the  combination 


(43) 


with 


p  E  (l  +  a''*'*)*”  -  a'*'' 


(44) 


From  (37),  (41),  and  (43),  the  expression  for  a  single  spectral  com¬ 
ponent  is 


F  =  A 


1/2  ^ 


(45) 


The  factor  A  is  chosen  to  match  Sudden's  result  for  z  «  when 
kQ/7r  is  large: 


.  iTr/4  1/2 
A  =  e  C 


dzq  ]  .  (46) 


Finally,  the  complete  expression  for  F  is 


-  Jj  <iSj  dSj  <;<s^, 


1/4 


(E^  -  1)»‘ 


X  e 


-kQ 


(l-ia-^^'>)o«p(-lU  /  JB*(f  .iVf  0 


X  exp  ^-ik(SjX  +  S2y)J  . 


(47) 
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This  is  the  desired  generalization,  holding  uniformly  throughout  an 
attenuating  layer,  of  Sudden's  result  for  a  single  reflection  level. 

ASYMPTOTIC  LIMITS 

The  expression  (47)  must  be  shovm  to  produce  the  correct  asymp¬ 
totic  limits  for  a  thick  evanescent  layer.  For  kQ/ir  large,  we  apply 
the  asymptotic  approximations  in  terms  of  Airy  functions  given  by 
Abrarnouitz  and  Sbegun  [1965]: 


(53) 
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From  (44) , 


1  -kQ 
y  cw  -  e 


(54) 


Combining  (43) ,  (48) ,  (49) ,  and  (54)  gives 

X  [2  Ai(-0  -  i  Y  Bi(-?)]  .  (55) 

For  z  <<  the  Airy  functions  can  be  replaced  by  their  asymptotic 
approximations ; 


Ai(-?) 


-1/2  ^-1/4  . 

TT  ?  Sin 


j  ■'T 

dzq  + 


(56) 


and 


Bi(-C) 


-1/2  ^-1/4 
Tr  f  cos 


(■/' 


j  ^  TT 

dzq  + 


(57) 


Expanding  the  trigonometric  functions  In  their  exponential  representa¬ 
tions,  substituting  into  (55),  and  collecting  terms,  we  obtain 

^  #  c)  “  (- !)  -  4^) 


X  exp 


("I  dz^  ^cxp  ^-ik  ^  d 

cp  ^-i2k  y*  dz^^l  -  —2—^  exp  ^ik  y*  dzqj[ 


(58) 
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Substituting  (58)  and  (46)  Into  (45)  gives 


This  represents  the  usual  decomposition  of  the  field  far  below  the 
layer  into  upgoing  and  downgoing  waves  with  reflection  coefficient 


(60) 


which  in  turn  yields 


I  „ 1 2  ,  -2kQ 
R  1  -  e  ^ 


1  +  e 


-2kQ 


(61) 


-kQ 

to  terms  of  second  order  in  e  .  The  last  result  is  a  generalization 
of  the  reflection  coefficient  expression  for  a  parabolic  layer  {Huihlcn, 
1961a]. 

For  z  near  the  higher  root  z^  of  q,  applying  (48)  and  (49)  to  (43) 

gives 


X  fBi(-C)  -  iAi(-C)]  ,  (62) 


NUMERICAL  KXAMPl.K 

The  energy  transmission  through  an  evanescent  layer  Is  given 
approximately  by  (61)  or  (67)  as 


(68) 
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For  the  hyperbolic  secant  layer  used  in  Sec.  IV  (witli  7.  =  0),  (36) 
gives 


Q  = 


i/C 


sech^u 


- 


C  <  C  , 
P 


=  0  ,  ^ 

where  C  =  a.  Rays  with  C  >  C  (or  S  <  S  )  are  refracted  but  not 
P  P  P 

attenuated  by  the  layer.  Following  Hudden  [1961/7],  the  parameter  a 

represents  f  /f,  where  f  and  f  are  the  wave  and  penetration  fre- 
P  P 

quencies,  respectively.  For 


k 


a 


2iif 


ca 


»  1 


a/C  1  and 


(70) 


Substituting  into  (68)  gives 

|T|2-exp[-2,to(f  -  l)]  ,  (71) 

Figure  5  plots  |t|  versus  S  computed  with  (71)  for  a  =  0.05  km  , 

k  =  4tt  km  and  various  values  of  a. 

2 

As  expected,  [t]  falls  rapidly  with  increasing  S  beyond  S^. 
Nevertheless,  for  frequencies  close  to  the  penetration  frequency 
(a  ~  1) ,  plane  wave  components  within  ~3  deg  of  vertical  will  suffer 
less  than  10  dB  loss  and  should  be  included  in  the  spectral  integral  (47). 
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VI.  CONCLUSIONS 

The  numerical  motliod  presented  here  for  modeling  radio  field 
strengths  in  caustic  and  cusp  reglon.s  lias  two  advantages  over  the 
asymptotic  method: 

1.  Use  of  a  ray  tracing  program  permits  phase  integrals  for 
complicated  media  to  be  computed  straightforwardly. 

2.  Use  of  numerical  qu.idrature  for  spectral  integration  obvi¬ 
ates  the  problems  of  locating  and  cah-ulating  the  contribu¬ 
tions  of  the  stationary  components  in  the  integral,  par¬ 
ticularly  when  higher  order  asymptotii-  methods  are  necessary. 

Numerical  Integration  is  both  feasible  and  much  easier  even  for  the 
simple  cases  discussed  here.  Our  algorithm  avoids  the  excessive  com¬ 
puter  running  time  often  associated  with  h.ighJy  oscillatory  integrands. 
For  simple  profiles  amenable  to  asymptotic  c.ilculat ion,  our  method 
yields  results  nearly  identical  to  those  of  Miit'iui  11976/;);  moreover, 
it  is  applicable  to  a  much  broader  r.ange  of  profiles.  Wi-ak  collisions 
can  be  included  through  use  of  attenuation  information  provided  by 
the  ray  trace.  Possible  general  I z;it ions  would  include  treatment  of 
anisotropic  media  and  leakage  through  thin  layers.  The  latter  prob¬ 
lem  is  partially  solved  by  extending  the  angular  spectral  integral 
to  include  regions  throughout  the  vicinity  of  a  thin  layer. 
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Appendix  A 

CURVE  FITTING  THE  RAY  INTERCEPT  DATA 


This  appendix  considers  the  problem  of  constructing  phase  inte¬ 
grals  from  ray  trace  intercept  data  on  a  plane  at  height  z.  Because 
of  the  polar  symmetry  in  S,  it  is  sufficient  to  consider  rays  in  the 
plane  y  =  0.  The  up-  and  downgoing  phase  integrals  z)  and 

z)  defined  by  (17)  and  (18),  respectively,  are  found  from  rays 
penetrating  the  plane  z  from  below  or  above.  Using  (16),  one  can  find 
<{i^  or  <fij  up  to  constants  by  integrating  an  analytic  representation  of 
the  ray  intercepts  x(S)  with  respect  to  S.  The  problem  then  is  one 
of  curve  fitting  the  ray  data  to  a  suitable  set  of  basis  functions 
over  the  range  of  S  values  that  produce  rays  intersecting  the  plane. 
For  the  upgoing  rays,  that  range  extends  from  zero  to  the  value 
for  which  the  ray  is  tangent  to  the  plane.  The  downgoing  range  runs 
up  to  from  S^,  which  labels  the  smallest  incidence  angle  ray  not 
penetrating  the  medium. 

For  a  given  height  z,  we  require  a  fit  to  both  branches  of  x(S) 
for  Sq  s  S  i  where  <  S^(z).  On  the  upgoing  branch,  s  0; 
on  the  downgoing  branch,  Sq  >  S^.  The  values  of  and  must 
be  chosen  to  include  all  rays  near  the  x  values  for  which  the  field 
is  to  be  computed.  Because  of  the  well-known  problem  of  numerical 
instability  in  least-squares  fitting  of  a  polynomial  of  degree  greater 
than  about  five,  orthogonal  polynomials  are  often  used  as  basis  func¬ 
tions.  A  particularly  convenient  choice  is  the  Chcbyshcv  polynomial 
set,  because  of  its  orthogonality  with  respect  to  summations  in 
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addition  to  its  ability  to  produce  a  uniform  fit  over  the  function 
domain.  Those  properties  are  well  described  by  Bahp.id  [1968].  To 
use  the  orthogonality  property,  rays  must  be  traced  at  definite,  dis¬ 
crete  values  of  S  given  by  the  equation 


I  ^  k  £  m  . 


(A. 6) 


For  m  =  N  -  1,  this  procedure  produces  o  collocation  polynomial  through 


the  data  points  x(S^);  for  smaller  m,  it  yields  a  least-squares  fit. 
Because  of  the  orthogonality  for  a  given  N,  the  coefficients  are 
fixed  for  all  fits  of  order  m  s  N  -  1.  Therefore,  in  contrast  to  the 
usual  least-squares  polynomial  fit,  increasing  the  order  of  the  approx¬ 
imation  does  not  result  in  a  reshuffling  of  lower  order  coefficients. 

We  illustrate  the  curve-fitting  procedure  using  data  produced 
from  the  tJonci'/!'iir;f>hf;ni',on  [1975]  ray  tracing  program  for  the  sech 
ionospheric  model  (30)  with  the  parameters  in  Sec.  IV  and  z  =  89  km. 

We  consider  licre  only  the  downgoing  rays  and  set  =  0.47,  Sj  =  0,63, 

N  =  12.  The  table  lists  the  coefficients  produced  for  each  fit  with 
0  ^  m  ^  11  and  the  associated  RMS  defined  as 


f  N-1 

P  m  - 

2' 

RMS  =  1 

x(Sj)  -23  Vk«>i> 

1  1=0 

L  k=o 

Only  the  highest  degree  is  given,  because  the  k  <  m,  are  already 
determined.  The  RMS  is  given  in  kilometers,  so  a  fifth-order  fit  pro¬ 
duces  a  representation  of  the  ray  intercepts  valid  to  about  0.03  km. 

The  method  thus  should  be  accurate  enough  to  produce  reliable  results 
in  the  angular  spectral  integration  when  used  with  Mdn/i.n'r,  frequency 
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COEKFICIENTS  PRODUCED  BV  CHEBYSHEV  FITS 
TO  RAY  TRACE  INTERCEPT  DATA 


Order 
of  Fit, 

m 

Highest  Degree 
Coefficient, 

(X 

in 

0 

91.143250 

J .8698 

1 

-2.358955 

0.8448 

2 

0.771790 

0.6449 

3 

0.901365 

0.0983 

4 

0.027617 

0.0963 

5 

0.128946 

0.031 1 

6 

-0.031670 

0.0218 

7 

0.015146 

0.0190 

8 

-0.002624 

0.0189 

9 

0.013373 

0.0165 

10 

-0.016300 

0.0143 

I J 

-0.019619 

0.0029 
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Appendlx  B 
PLP  ALGORITHM 


This  appendix  describes  the  numerical  quadrature  algorithm  of 
Woodze  [1976]  and  Bavakat  [1976]  used  in  performing  the  plane-wave 
angular  spectral  integrations.  The  spectral  Integrals  are  of  the  form 


l(x) 


S 


S  . 


min 


dS  g(S) 


^-ikP(S,x) 


(B.l) 


where  both  g(S)  and  P(S,  x)  are  continuous  functions  over  the  inter¬ 
val  S  .  <  S  <  S  .  Because  k  can  be  large,  the  Integrand  contains 

min  max 

a  highly  oscillatory  phase  term.  Since  P(S,  x)  is  continuous,  we  can 
subdivide  the  integration  Interval  into  N  segments  such  that  in  the 
ith  segment 

(S  -  S.)  ,  (B.2) 


P(S,  x) 


P(S  x)  +  — 


where  S  =  S  ,  +  AS(i  -  1/2),  i  =  1,  N  +  1,  and  AS  =  (S  -  S  .  )/N. 

i  min  mm 

It  is  also  assumed  that  N  is  large  enough  that  g(S)  =“  g(Sj)  within 
each  element.  Substituting  into  (B.l)  produces 


I(x) 


N+1 


-lkP(Sj,x) 

e 


S.+AS/2 

/ 

S^-AS/2 


-ik(:(s  x)(.s-s  ) 

dS  e 


N+1 


AS 


s(s^) 


i  =  l 


-ikP(Sj,x)  .sin  kG(Sj,  x)  AS/2 
'  kGTs'jV^  AS72  ~  ’ 


(B.3) 


-aV.,... 
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where  (:(S^,  x)  =  g— j  .  Because  of  the  approximation  (B.2),  this  pro- 
cedure  has  been  named  the  piecewise  linear  phase  (PLP)  algorithm. 

Its  extension  to  two  dimensions  is  straightforward. 
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SYMBOLS 


A(S,  z)  •  Budden's  reflection  region  factor,  given  by  (4)  or 
its  asymptotic  limit  (5). 

Ai(5),  Bi(5)  =  Airy  integral  functions. 

C  =  z-direction  cosine  of  plane  wave  component  in  free 
space . 

* 

E  (a,  x)  =  parabolic  cylinder  function. 

F  =  component  of  electromagnetic  field. 
g^(S)  =  nt,}/  azimuthal  «-omponent  ol  (I. 

G(S^,  S^)  =  spectral  weighting  function. 

H  =  Hamiltonian. 

J  =  nth  Bessel  function. 

n 

k  =  free-space  wave  number,  u/c,  where  o)  Is  the  angular 
frequency  and  c  is  the  speed  of  light. 
n(z)  =  refractive  index. 

q  =  [n  (z)  -  S  ] 

Q  =  integral  of  (q|  through  the  evanescent  region. 

S  =  (S^  + 

S^,  S2  =  X-  and  y-direction  cosines  of  plane  wave  component  in 
free  space. 

W(a,  x)  =  Weber's  function. 

X,  y,  z  =  Cartesian  coordinates  with  the  z-axis  vertically 
upwards . 

Z  =  v/(i),  where  V  is  the  electron  collision  frequency  and 
0)  Is  the  angular  frequency. 
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5  =  argument  of  Airy  function  given  by  (3). 


,  2  .  2.1/2 
(x  +  y  ) 


T  =  independent  variable  in  Hamilton’s  ray  equations. 
<})  =  tan  ^(82/3^) . 

^2’  ~  pluise  integral  for  upgoing  wave  at  height  z. 

^2’  se  integral  for  downgoing  wave  at  height  z. 

tp  =  tan  ^  (y/x) . 
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